#1. Packages ----
library(ncdf4)
library(ncdf4.helpers)
library(PCICt)
library(sf) 
library(tidyverse)

library(ggplot2)
library(ggmap)
library(usmap)
library(mapdata)
library(ggthemes)
library(dplyr::select())
library(RColorBrewer)
library(viridis)
library(maps)
library(cowplot) # for plot_grid() function - good to arrange multiple plots into a grid
library(grid)
library(gridGraphics)
library(patchwork)
library(gridExtra)
library(colorspace)
#library(raster)
library(rgdal)
library(broom)
library(exactextractr)

library(wesanderson)
library(MetBrewer)
library(sp)
library(sf)
library(classInt)


library(ggpubr)
require(maps)
require(viridis)

theme_set(
  theme_void()
)

theme_set(theme_classic())


#2.Paths ----


BASEPATH="C://Users//wb520443//Dropbox//Research Projects//Global pollution//2_analysis//Replication"


#3. Data import ----
world <- read_rds(paste0(BASEPATH,"/figdat/world_borders_WB.rds"))



global_pm_income <- read_csv(paste0(BASEPATH,"/figdat/global_pm_income.csv"))



sum(global_pm_income$ISO %in% world$ISO3) #166 countries match
sum(global_pm_income$ISO %in% world$ISO3 == F) #90 do not
global_pm_income$ISO[global_pm_income$ISO %in% world$ISO3 ==F] #print countries that don't match (all small islands)

plotdat <- global_pm_income %>% select(ISO3 = ISO, starts_with("mean")) 

world@data$plotOrder <- 1:nrow(world@data)

world@data$ISO3<-world@data$ISO3166_1_

world@data <- left_join(world@data, plotdat) %>% arrange(plotOrder)



#4. PM Map ----


pal <- met.brewer("OKeeffe2") 
int <- classIntervals(world@data$mean_pm, style = "fixed", fixedBreaks = c(0, seq(0, 50, 6), 100))
col <- findColours(int, pal)



pdf(paste0(BASEPATH,"//figures//map-global-pm.pdf"), width = 7.2, height = 4)
par(mar = c(0,0,0,0))

par(bg = 'white')
plot(world, col = 'gray90', border = NA, ylim = c(-85, 75), xlim = c(-170, 170))
plot(world, border = NA, col = col, add =T)
text(x = coordinates(world)[,1], y = coordinates(world)[,2], label = world@data$ISO3, cex = 0.2, col = 'gray90')



plotrix::gradient.rect(-115,-85,115,-72.5,
                       col=plotrix::smoothColors(colorRampPalette(pal)(10)), 
                       gradient="x", border = 'black')

text(x = seq(-105,105,210/9) ,y = -90, labels = c(5,10,15,20,25,30,35,40,45,50), cex = 1.1)
segments(x0 = seq(-118.1,-117.8,.03),y1 = 33.465, y0=33.46)
text(x = 0, y = -65, labels = "Average annual PM2.5",cex=1.3)


dev.off()


#5. Per cap GDP map ----



pal <- met.brewer("Hokusai2") 
int <- classIntervals(world@data$mean_gdp_per_cap, style = "fixed", fixedBreaks = c(0, 2000, 4000, 6000, 8000, 10000, 12000, 14000, 16000, 18000))
col <- findColours(int, pal)



pdf(paste0(BASEPATH,"//figures//map-global-income-pm.pdf"), width = 7.2, height = 4)
par(mar = c(0,0,0,0))

par(bg = 'white')
plot(world, col = 'gray90', border = NA, ylim = c(-85, 75), xlim = c(-170, 170))
plot(world, border = NA, col = col, add =T)
text(x = coordinates(world)[,1], y = coordinates(world)[,2], label = world@data$ISO3, cex = 0.2, col = 'gray90')



plotrix::gradient.rect(-115,-85,115,-72.5,
                       col=plotrix::smoothColors(colorRampPalette(pal)(9)), 
                       gradient="x", border = 'black')

text(x = seq(-105,105,210/8) ,y = -90, labels = c(2, 4, 6, 8, 10, 12, 14, 16, 18), cex = 1.1)
segments(x0 = seq(-118.1,-117.8,.03),y1 = 33.465, y0=33.46)
text(x = 0, y = -65, labels = "Average per capita income, 2016-2020, 2015USD, 000s",cex=1.3)


dev.off()


#6. Urban map ----



# pal <- met.brewer("Degas") 
# int <- classIntervals(world@data$mean_pct_urban, style = "fixed", fixedBreaks = c(0, 12, 24, 36, 48, 60, 72, 84, 96))
# col <- findColours(int, pal)

pal <- colorRampPalette(col =c("#E3E6B8","#EFE085","#F5BD6C","#D37956","#910F3C"))(256)
int <- classIntervals(world@data$mean_pct_urban, style = "fixed", fixedBreaks = seq(10,100,10))
col <- findColours(int, pal)




pdf(paste0(BASEPATH,"//figures//map-global-urban.pdf"), width = 7.2, height = 4)


par(mar = c(0,0,0,0))

par(bg = 'white')
plot(world, col = 'gray90', border = NA, ylim = c(-85, 75), xlim = c(-170, 170))
plot(world, border = NA, col = col, add =T)
text(x = coordinates(world)[,1], y = coordinates(world)[,2], label = world@data$ISO3, cex = 0.2, col = 'gray90')



plotrix::gradient.rect(-115,-85,115,-72.5,
                       col=plotrix::smoothColors(colorRampPalette(pal)(10)), 
                       gradient="x", border = 'black')

text(x = seq(-105,105,210/9) ,y = -90, labels = seq(10,100,10), cex = 1.1)
segments(x0 = seq(-118.1,-117.8,.03),y1 = 33.465, y0=33.46)
text(x = 0, y = -65, labels = "Percent of the population living in urban setting",cex=1.3)


dev.off()

#7. Scatter plot ----


scatter<-ggplot(global_pm_income, aes(y = mean_pm, x=mean_gdp_per_cap )) +
  geom_point(size=2, colour = "grey70") +
  geom_text(aes(label=ifelse(ISO=="QAT",as.character(ISO),'')),hjust=1.12,vjust=1.1) +
  geom_text(aes(label=ifelse(ISO=="ARE",as.character(ISO),'')),hjust=-0.1,vjust=1.1) +
  geom_text(aes(label=ifelse(ISO=="MAC",as.character(ISO),'')),hjust=1.1,vjust=1.1) +
  geom_text(aes(label=ifelse(ISO=="MCO",as.character(ISO),'')),hjust=1.1,vjust=0.1) + 
  geom_smooth(method = "loess", size = 1.5, color='blue',level=0.0) +
  labs(color = "Percent of country\nthat is urban")+
  theme(legend.direction="vertical",legend.justification=c(1,1), legend.position=c(1,1))+
  theme(text = element_text(size = 15))+ 
  theme(legend.background=element_blank())+
  theme(axis.line.x = element_blank(),
        axis.line.y = element_blank())+
  xlab("Average per capita income,\n2016-2020, 2015USD")+
  ylab("Average PM2.5")+
  labs(subtitle = "") +
  theme(plot.subtitle = element_text(size = 30,  hjust = 0.02))



ggsave(filename = paste0(BASEPATH,"//figures//mean_pm25_vs_income_global.png"), 
       plot = scatter,
       width = 5,
       height = 6,
       units = c("in"),
       dpi = 300)


